NPSOR-91-030 

NAVAL  POSTGRADUATE  SCHOOL 

Montepej,  California 


* 


I 


A  KALMAN  FILTER  FOR  A  POISSON  SERIES 
WITH  COVARIATES  AND  LAPLACE 
APPROXIMATION  INTEGRATION 


Donald  P.  Gaver 

// 

Patricia  A.  Jacobs 


September  1991 


Approved  for  public  release;  distribution  is  unlimited. 
Prepared  for: 

Naval  Postgraduate  School 
Monterey,  CA  93943 


r 

FedDocs 
D  208.14/2 

NPS-OR-9 1-030 


I  fr( 

Loo^  /4(^' 

-  op- OdO 


f'Af  T 

SCHtXflu 

.  4  J-,a4S-x»0»JPi» 


NAVAL  POSTGRADUATE  SCHOOL, 
MONTEREY,  CALIFORNIA 


Rear  Admiral  R.  W.  West,  Jr.  Harrison  Shull 

Superintendent  Provost 


This  report  was  prepared  in  conjunction  with  research  funded  under  the  Naval 
Postgraduate  School  Research  Council  Research  Program. 


This  report  was  prepared  by: 


UNCLASSIFIED 


.  ,rY  X  ...K. 

•AYAL  POSTGRADUATE  SCHO<^- 

CAXiIFORNtav 


SECURITY  CLASSIFICATION  OF  THIS  PAGE . 


DUDLEY  KNOX  LIBRARY 
NAVAL  POSTGRADUATE  SCHOOL 
MONTEREY  CA  93943-5101 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
OMBNo.  0704-0188 


la.  REPORT  SECURITY  CLASSIFICATION 

UNCLASSIFIED 


1b.  RESTRICTIVE  MARKINGS 


2a.  SECURITY  CLASSIFICATION  AUTHORITY 


2b.  DECLASSIFICATION/DOWNGRADING  SCHEDULE 


3.  DISTRIBUTION  /AVAILABILITY  OF  REPORT 

Approved  for  public  release;  distribution  is 
unlimited. 


4.  PERFORMING  ORGANIZATION  REPORT  NUMBER(S) 

NPSOR-91-030 


5.  MONITORING  ORGANIZATION  REPORT  NUMBER(S) 


6a.  NAME  OF  PERFORMING  ORGANIZATION 

Naval  Postgraduate  School 


6b.  OFFICE  SYMBOL 

(If  applicable 

OR 


7a.  NAME  OF  MONITORING  ORGANIZATION 

O&MN,  Direct  Funding 


6c.  ADDRESS  ( City,  State,  and  ZIP  Cod^ 

Monterey,  CA  93943 


7b.  ADDRESS  (Gty.  State,  and  ZIP  Codd) 

Monterey,  CA  93943-5006 


8a.  NAME  OF  FUNDING/SPONSORING 
ORGANIZATION 

Naval  Postgraduate  School 


8b.  OFFICE  SYMBOL 

(If  applicable) 


9.  PROCUREMENT  INSTRUMENT  IDENTIFICATION  NUMBER 


8c.  ADDRESS  ( City,  State,  and  ZIP  Code 

Monterey,  CA  93943 


10.  SOURCE  OF  FUNDING  NUMBERS 


PROGRAM 

PROJECT 

TASK 

WORK  UNIT 

ELEMENT  NO. 

NO. 

NO. 

ACCESSION  NO. 

11.  TITLE  (Include  Security  Classification) 

A  Kalman  Filter  for  a  Poisson  Series  with  Covariates  and  Laplace  Approximation  Integration 


12.  PERSONAL  AUTHOR(S) 

Donald  P.  Gaver  and  Patricia  A.  Jacobs 


13a.  TYPE  OF  REPORT 

Technical  Report 


13b.  TIME  COVERED 
FROM TO 


14.  DATE  OF  REPORT  (Year,  month  day) 

September,  1991 


IS.  PAGE  COUNT 

47 


16.  SUPPLEMENTARY  NOTATION 


17.  COSATI  CODES 

FIELD 

GROUP 

SUB-GROUP 

1 8.  SUBJECT  TERMS  (Continue  on  reverse  if  necessary  and  identify  by  block  number) 

Poisson  time  series  with  covariates;  hierarchical  model;  Kalman 
filter;  Laplace  method 


1 9 .  ABSTRACT  (Continue  on  reverse  if  necessary  and  identify  by  block  number) 

A  hierarchical  model  for  a  Poisson  time  series  is  introduced.  The  model  allows  the  mean  or 
rate  of  the  Poisson  variables  to  vary  slowly  in  time;  it  is  modeled  as  the  exponential  of  an  AR/1 
process.  In  addition  the  rate  is  influenced  by  a  covariate.  The  Laplace  method  is  used  to 
recursively  update  some  model  parameter  estimates.  Frankly  heuristic  methods  are  explored  to 
estimate  other  of  the  underlying  parameters.  The  methodology  is  checked  against  simulated 
data  with  encouraging  results. 


20.  DISTRIBUTION/AVAILABILITY  OF  ABSTRACT 

0  UNCLASSIFIED/UNLIMITED  QSAMEASRPT.  Q  DTIC  USERS 

2 1 .  ABSTRACT  SECURITY  CLASSICIATION 

UNCLASSIFIED 

22a.  NAME  OF  RESPONSIBLE  INDIVIDUAL 

D.  P.  Gaver 

22b.  TELEPHONE  (/nc/ud©  Area  Code; 

(408)  646-2605 

2c.  OFFICE  SYMBOL 

OR/Gv 

DD  Form  1473,  JUN  86 


Previous  editions  are  obsolete. 


SECURITY  CLASSIFICATION  OF  THIS  PAGE 


I 


I 


4 


I 

t 


I 


j 


A  KALMAN  FILTER  FOR  A  POISSON  SERIES  WITH  COVARIATES 

AND  LAPLACE-APPROXIMATION  INTEGRATION 

D.  P.  Gaver 
P.  A.  Jacobs 


0.  INTRODUCTION 

The  Poisson  model  is  an  initial  idealized,  but  plausible,  off-the-shelf  tool 
for  representing  point-process  data  of  nearly  any  kind,  cf.  Feller  (1966)  and 
Cox  and  Lewis  (1968).  However,  to  be  more  descriptive,  and  even  predictive, 
representation  of  non-homogeneity  in  space  or  time  may  be  needed.  For 
instance  the  occurrence  of  rare  events  in  space,  such  as  the  occurrence  of 
extreme  heights,  i.e.  above  a  level  of  Arctic  ice  along  a  transsect  or  encounters 
with  ice  keels  along  a  submarine  track  at  constant  (deep)  depths,  may  well 
appear  roughly  Poisson.  A  better  description  of  these  events  may  require 
more  detail  than  a  simple  mean  or  rate:  some  account  of  regional  and 
seasonal  variation  could  be  needed  for  true  accuracy;  for  instance  the 
intervention  of  natural  gaps  along  a  path  in  Arctic  ice  will  occur  if  the  latter 
crosses  leads  (open  water  amidst  an  ice  pack).  For  another  example,  demands 
for  spare  parts  in  a  logistics  system  often  are  roughly  Poisson,  or  compound 
Poisson,  but  with  a  mean  or  rate  that  changes  erratically  but  slowly  in  time. 
Demands  for  communication  and  computer  facilities  exhibit  a  similar 
temporal  pattern,  and  there  are  a  great  many  other  examples.  To  summarize, 
variation  in  a  fundamental  Poisson  rate  or  mean  is  very  often  encountered  in 
practice;  it  is  possible  that  if  this  variation  is  slow-moving  or  persistent 
enough  it  can  be  exploited  for  short-term  forecasting. 
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In  this  paper  we  study  a  model-based  procedure  for  forecasting  in  the 
kind  of  environments  described.  The  model  introduced  allows  the  mean  or 
rate  of  the  Poisson  process  to  be  itself  a  random  process;  the  exponential  of  an 
AR/1  autoregressive  process.  In  addition  the  rate  is  influenced  by  a 
covariate.  We  then  recursively  update  the  parameter  estimates  using  an 
approximation  based  on  the  Laplace  method;  cf.  de  Bruijn  (1958).  The  approach 
resembles  that  of  Delampady,  Yee  and  Zidek  (1991),  but  frankly  heuristic 
methods  are  used  to  estimate  certain  of  the  underlying  parameters.  The 
methodology  is  checked  against  simulated  data  with  encouraging  results. 

1.  FORMULATION 

Consider  the  following  Poissonian  model  for  count  data: 

|!" 

Vf 

where 

lif  =  “Ilf -1  +  (1-2) 

with  {©j}  independent  normal/Gaussian  random  variables  with  mean  0  and 
variance  IN ^  independent  of  {p,;;  i  <  f-1},  {Y,-;  i  <  t-l],  {xj,  {h^}  and  p.  Another 
hierarchical  time  series  model  for  count  data  can  be  found  in  Harvey  et  al. 
(1989). 

The  purpose  of  this  paper  is  to  suggest  a  Kalman  filter-like  procedure  to 
produce  successive  estimates  of  P  and  Pj  as  new  data  becomes  available.  The 
procedure  is  based  on  a  Laplace  approximation  to  an  integral.  A  Bayesian 
approach  to  a  similar  problem  is  being  investigated  by  Delampaday  et  al. 
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(1991).  A  Bayesian  approach  to  time  series  can  be  found  in  West  et  al.  (1989) 
and  for  a  more  recent  computational  approach  in  Carlin  et  al.  (1991). 


2.  AN  APPROXIMATE  UPDATING  PROCEDURE 

Assume  that  the  posterior  distribution  of  (P,  given  {y„  i<  t-1}  is 

bivariate  normal  with  mean  (b,_j,  wi,_j),  Var[p]  =  Var[Hj_J  =  C^_^  and 

Corr[p,n,_i]  = 

Since  it  is  known  that 

\JLt  =a\Lt.i+<Ot, 


the  prior  distribution  of  (P,pO  is  bivariate  normal  with  mean  octn^_^), 
Var[p]  =  Var[n,]  =  R,  =  oc^C^_^  +  W,  and  Corr[p,n,]  =  r,  =  ap,_j  yjQ-i/Rt. 

The  forecast/prediction  distribution  of  Y;  in  terms  of  data  up  to  t-1  and 
the  covariate  value  at  t  is 


P  {Y,  =  i/,|Y,-  =  yj  J<t-  <  t,  hi  ,i  <  f} 


=  j*  expj- 


exp 


Vt'-  2n^Tt.iRt[l  -  r^) 

(fc  -  fcf )  2rt[b  -  bf)(z  ■  wf)  (z  -  mf) 


Tf-i 


dzdb.  (21) 


where  bf*  =  bf.i  and  mf  =  anti  .j. 


We  now  approximate  the  integral  by  the  Laplace  method;  cf.  Easton 
(1991),  Cox  and  Hinkley  (1974),  de  Bruijn  (1958).  Let  the  exponent  of  the 
integrand  be 
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(2.2) 


g{b,z)  =  +  yt[Xtb  +  z]  +  ytlnht  -  ^ 


I,.i  V?Ta'  "  R, 


+  K 


where  J<C  is  a  constant. 


±S(b,z)  -  y,  -  e>‘''*^h,  -  ^ 


02 


1-r/ 


K'^f)  .  (^-^D 
^  -1 


=  (y,  -  -  [l  - 


-1 


(^  •  •  o 

■^f-i  V  - 1 


afc‘ 


-^(fc,2)  =  -[i 


■1  1 


^f-i 


awz 


^(b,2)  =  +{r^)[l  ■  r^)  ^Tj.iR,) 


05 


Use  a  Newton  procedure  to  solve  the  system  of  equations 

-1^ 


O-A, 

db^ 


for  mt  and  bt-  Solve  for  T(,  Q  and  pt  using  the  relation 


(2.3) 

(2.4) 

(2.5) 

(2.6) 

(2.7) 

(2.8) 
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'Tf 

r 

db^ 

mt,bt 

dbdz 

mM 

\ - 

Ct 

1 

,dbdz  ^ 

mt.bt 

& 

mt,bt  ^ 

The  posterior  distribution  for  (P,p/)  given  yt,  Dt-i  is  approximated  by  a 
bivariate  normal  distribution  with  mean  (bt,  rrit)  and  variance-covariance 
matrix 


(2.10) 


2.1  Summary  of  the  Newton  Procedure 

0.  Start  with  estimates  of  the  parameters  of  the  prior  bivariate  distribution 
of  (P,Pf);  that  is,  estimates  of  the  mean  (bt-i,  anit-i),  Corr(P,Pf)  =  r/,  Var(P)  =  T/_i 
Var[p/]  =  Rf.  Set  bf  =  b^^,  wi*f  =  am^_y 

1.  Solve  the  sy stern  of  linear  equations  (in  b  and  m) 


for  b,  m . 

2.  If maxj^(|b®  -b  -m|/m)j<  0.001,  set  mt  =  m  and  bt=  b  and  go  to 

3.  Otherwise  set  b^  =  b  and  m^t  =m  and  return  to  Step  1  unless  Step  1  has  been 
returned  to  49  times  in  which  case  set  nit  =  m^t  arid  bt  =  bf;  go  to  3. 

3.  Return  (bt,  nit)  as  the  estimate  of  the  posterior  mean  of  the  bivariate 
normal  distribution  of  0,Pf). 
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4.  To  find  the  variance-covariance  matrix  of  the  posterior  distribution  of 
(P,pt)  evaluate 


db‘ 


idzdb 


dzdb 

:.2 


g{h>^t) 


■^gibt>^) 

dz 


set  it  equal  to 

and  solve  for  t<  ,  pt  and  C^ 

2.2  Summary  of  Kalman  Procedure 

In  summary  the  approximate  Kalman  procedure  is  as  follows 
0.  Start  with  the  parameters  of  the  (approximate)  posterior  bivariate 
normal  distribution  of  0/lit-i)  having  mean  (bt-\,  Corr(p,p,(_i)  =  pt-\, 


Var[p]  =  Tt_i  Var[^,_i]  =  Cm- 


1.  Update  Corr(p,p.f_i)  and  Var[p,(.i]  using  (1.2)  to  obtain  the  prior 
bivariate  normal  distribution  of  (p,|it)  having  mean  (^M/  otm/.i),  Var[p]  =  Tm/ 
Var[[it]  =Rt  =  a^Ct-\  +  Wt,  Corr(p,Hi)  =  r,  =  apt-i  >/q77^- 


2.  Observe  the  Poisson  count  x/i. 

3.  Invoke  the  Newton  procedure  of  Section  2.1  to  obtain  estimates  of  the 
parameters  of  the  approximate  posterior  distribution  of  (p,lit)  given  past 
observations  and  the  new  observation 


6 


To  obtain  moments  for  the  Poisson  mean  X(  =  fi<exp{XiP+^(}  note  that  since 
the  distribution  of  (fit,  pO  is  being  approximated  by  a  bivariate  normal 
distribution,  the  posterior  moments  of  Xt  are  approximately 

E[Xf  1  »  exp  xtbt  +mt+^[xjrt  +Ct +2xtptJ^J^\^ht 

Var[Xf  ]  “  hf  exp|2(3C^fcj  +  ntt)  +  [x(  +  Q  +  Ix^pt-J^ 

^[exp[xh  +Ct+2XtPtJ^J^}  -l\. 

The  estimate  of  Xfis  Xf  = 


2.3  A  Simulation  Example 

In  the  example  P  =  0.5,  {x^}  takes  the  values  {0.25,  .5, 1}  over  and  over;  ht=i. 

{cO/}  are  iid  normal  with  mean  0  and  variance  0.25  and  a  =  0.5.  Given  p  and  |i/, 

Yf  has  a  Poisson  distribution  with  mean  Xt  =  .  The  simulation  starts 

with  po  drawn  from  the  stationary  distribution  of  {p.,}  a  normal  distribution 

2  1 

with  mean  0  and  variance  W/(l-a  )  =  g.  Successive  p/are  computed  as 

Ilf  =  aiif-l  +®f; 

the  Xt  is  computed  and  a  Poisson  random  number  with  mean  X/  is  then 
generated.  The  simulation  data  are  the  Poisson  counts  and  (Xf);  f  =  1, ...,  100. 
The  random  numbers  were  generated  using  LLRANDOM  II;  cf.  Lewis  et  al. 
(1981). 

At  time  0,  the  Newton  procedure  is  initialized  at  po  =  0,  P  =  0,  po  =  ro  =  0, 
a  =  0.5,  W  =  0.25,  Cq  =  0,  rj  =  0.25,  tq  =  1;  note  the  a  and  W  are  assumed 
known.  The  data  point  yi  is  observed  and  the  Newton  procedure  is  used  to 
find  the  posterior  moments  [wii,  b  i,  pi,  Q,  Ti]. 
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The  da  a  point  1/2  is  the  observed  and  the  Newton  procedure  is  started 
with  [mi,  &  1,  ri.  Cl  +  Ti]  and  used  to  find  [m2,  b2,  pi,  C2,  T2],  etc. 

Results  of  the  simulation  appear  in  Figures  1-3.  In  Figure  1  the  count  data 
y^  appears  along  with  the  true  (dotted  line)  and  the  estimated  Aj  (solid  line). 
In  Figure  2  the  true  value  of  }i^  (dotted  line)  and  estimated  value  of  (solid 
line)  appear.  The  estimated  values  of  the  standard  deviation  of  fi^,  ^/c,  ,  also 
appear  in  Figure  2.  The  estimated  value  of  p's  estimated  standard 
deviation  yjr^,  and  the  estimated  value  of  appear  in  Figure  3.  Figures  1-3 
indicate  that  the  procedure  performs  satisfactorily.  The  apparent  oscillation 
in  some  of  the  figures,  (particularly  is  due  to  the  cyclic  nature  of  [jcJ.  Not 
surprisingly,  the  estimated  values  of  Xj  and  Pj  are  less  variable  than  the  true 
values.  They  appear  to  be  practically  acceptable. 

3.  APPROXIMATE  KALMAN  PROCEDURES  WHICH  INCORPORATE 
ESTIMATES  OF  THE  AUTOREGRESSIVE  PARAMETERS  a  AND  W: 
NAIVE  MOMENT  ESTIMATORS 

In  this  section  we  will  assume  that  {cOj}  are  independent  identically 
distributed  normal/Gaussian  random  variables  with  mean  0  and  variance  W 
with  cOj  independent  of  (jl-;  i  <  t-1],  (Y-;  i<  t-  1},  [h^]  and  p. 

3.1  Estimates  of  the  Autoregressive  Parameters  a  and  W 

If  {p,;  t  <  T}  were  observable,  then  one  could  estimate  a  and  W  by 
maximum  likelihood;  that  is,  the  likelihood  function  is 


(3.1) 


or 
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(3.2) 


T 

lnL(W ,a)  ^  liW ,a)  = 

^=1  ^ 

T  T 

— InVV  — ln2;r. 

2  2 


Now 


gives 


and 


from  which 


da 


T 

■Ife 


t=l 


=  0 


T  T 

a(T)  =  - 1  /  Z 


f=i  f=i 


a/ 

aw 


1  1 
2w2 


Z(Mf  -  m-if 

f=l 


1^=0 

2  W 


(3.3) 


(3.4) 


(3.5) 


(3.6) 


Unfortunately  is  not  observable.  One  possible  estimate  of  is  the 
posterior  mean  m^  of  subsections  (2.1)  and  (2.2).  In  this  case  the  corresponding 
estimators  of  a  and  W  are 


^miT)  = 

f=l  f=l 


(3.7) 


W„(T)-iXK  ■d„(T)m,.i)^ 


(3.8) 


f  =  l 
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Another  estimate  of  ji^is 


My)  =  H  yt  + 


-  Xfbt 


(3.9) 


where  is  the  mean  of  the  approximate  posterior  distribution  of  P  given 
{y^;  s  <  f}.  The  corresponding  estimates  of  a  and  W  are 

T  T 

=  KMf(y)M/-i(y))/EM/(y)^  (3.io) 

f=i  i~\ 

and 

T 

Wy(T)  -  [My)  ■  <iy(3-)/x,  .i(y))  ■  (311) 

The  estimates  (3.7)  -  (3.11)  can  be  recomputed  at  every  time  T.  Other 
similar  estimates  can  be  obtained  by  not  recomputing  the  estimates  at  every 
time  T;  for  example,  choose  an  integer  5  >  1,  put 

a(t;5)  =  a((n  - 1)5;5)  (3.12) 

and 

VV(f;5)  =VV((n  -  1)5;5)  (3.13) 

if  (w  -1)5  <  t  <  n5; 

nS  n5 

a(M5;5)  =  (3.14) 

f=l  f=l 


1  >i5  2 

IV(m5;5)  =  a(«5;5)/i,.i)  (3.15) 


where  (1^  is  an  estimate  of  Pj.  Another  possibility  is  to  use  a  window  of  times 
to  compute  estimates  of  a  and  W. 
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The  following  is  a  summary  of  the  Kalman  procedure  of  Section  2  with 
the  addition  of  estimation  of  a  and  W. 

3.2  Summary  of  the  Kalman  Procedure  with  Estimation  of  a  and  W 

0.  Start  with  the  parameters  of  the  (approximate)  posterior  bivariate 
normal  distribution  of  (P,  having  mean  Corr(P,ji,_j)  =  p^_ 

\f  Var[p]  —  Var[}rj_j]  — 

1.  Update  Corr(P,  and  Var[|ij_j]  using  (1.2)  and  the  current  estimates 

of  a  and  W  to  obtain  the  prior  bivariate  normal  distribution  of  (P, 

2 

having  mean  {b^_y  Var[p]  =  Var[pJ  ^R^=  +  VJ ^_y 

Corr(P,p,)= 

2.  Observe  the  Poisson  count  y^. 

3.  Invoke  the  Newton  procedure  of  Subsection  2.1  to  obtain  estimates  of 
the  parameters  of  the  posterior  bivariate  normal  distribution  of  (P,  Pj) 

given  past  observations  and  the  new  observation  y^. 

4.  Compute  new  estimates  of  a  and  W: 

5.  Return  to  1. 

3.3  Results  of  Simulation  Experiments 

Various  simulation  experiments  were  carried  out  using  the  procedures 
outlined  in  this  section  to  estimate  a  and  W.  One  striking  phenomenon  was 
that  if  the  count  data  has  a  long  run  of  zeroes,  then  the  Kalman  procedure 
reacts  by  estimating  a  >  1  and  trying  to  estimate  Pj  to  be  a  negative  number 
large  in  absolute  value.  This  tendency  caused  the  numerical  Newton 
procedure  discussed  in  Section  2  to  become  unstable,  and  rendered  the 
predicted  value  of  very  slow  to  respond  to  positive  counts  when  they  did 
eventually  occur.  This  behavior  was  ameliorated  by  arbitrarily  putting  a 
lower  bound  of  -2  on  fi^.  A  run  of  zero  counts  could  also  make  the  estimate  of 
W  become  close  to  zero.  A  small  value  of  W  makes  the  Kalman  filter 
unresponsive  to  count  changes.  This  behavior  was  mitigated  by  recomputing 
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estimates  of  a  and  W  every  5  =  10  time  units  rather  than  at  every  time,  and 
using  all  previous  times  t  as  in  (3.7)-(3.8).  The  estimates  of  a  and  W  were  also 
not  computed  if  the  last  10  observed  counts  were  zero.  Further  a  lower 
bound  of  0.1  was  arbitrarily  set  on  estimates  of  W.  Such  heuristics  are  not 
claimed  to  be  optimal,  but  are  needed  to  allow  the  procedures  to  behave 
suitably.  A  search  for  more  systematic  procedures  can  be  carried  out  in 
future  work. 

The  following  are  empirical  observations.  The  estimates  of  a  and  W  using 
//(f)  =  nij  tends  to  produce  practically  acceptable  estimates  of  a  but 
underestimates  W — it  tends  to  make  it  very  small.  This  behavior  is  not 
unexpected  since  m^can  be  thought  of  as  a  smoothed  estimate  of  the  p.,  and  so 
will  have  a  smaller  variance  than  \L^.  This  underestimation  of  W  is  not 

.  Neither  does  using  "windows"  of  length  5 


improved  by  using  //^  =  In 


y/+2 


seem  to  solve  the  problem.  Other  procedures  for  estimating  a  and  W  are 
explored  in  the  following  two  sections. 


4.  ESTIMATION  OF  a  AND  W  FOR  THE  POISSON/NORMAL  MODEL: 
THE  FREEMAN-TUKEY  TRANSFORMATION 

This  section  reports  another  possible  approach  to  the  estimation  of  the 
autoregressive  parameters  for  a  Poisson/normal  model. 

The  model  is  as  before 


lif+1  =  m  +a)f+i 


(4.1) 


where  {ca^}  are  iid  normal  with  mean  0  and  variance  W; 

•  yUi  ,^i]  •  exp[-  - - -  .  (4.2) 

y  • 
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For  simplicity  of  notation  we  will  assume  h^=  1. 

If  W  and  a  are  known,  then  an  approximate  Kalman  procedure  has  been 
given  previously.  The  issue  here  is  the  estimation  of  a  and  W. 

Let 

^  g{x)=y/x  +  y/x  +  1  (4.3) 

and 

Z,  =  g(Y,).  (4.4) 

The  conditional  distribution  of  given  is  approximately  normal  with 

r  X 

mean  k{\L^  +  x^P),  where  /c(x)  =  [4e  +  Ij  ,  and  variance  1;  c.f.  Freeman  and 
Tukey  (1949,  1950)  and  Bishop,  Fienberg  and  Holland  (1975).  Thus  if  are 
iid  standard  normal,  given  =  (Yq,  Yj,  ...,  Y^) 
d 

Z(+i  +1  ^{+iP) 

=  k{a]ii+  (0{+i  +  XfViP)  +  Vf +1 

=  k{amt  +  a{\Lt  -  nit)  +  -  ^f))  + 

»  k{amt  +  Xt^ibt) +k'[amt  +  Xt+ibt)[a{\it  -  m^)  +x,+i(p  -  b^)  +  (D(+i]  +  V,+i.  (4.5) 
Hence, 

E[Zi+i|Df]  «  k{amt  +  Xt^ibt)  (4.6) 

Var[Zf+i|Df]  =  [fc'(am,  +  X(+ifc,)p[a^C,  +  Xf^+iT,  +  2aX(+iPf  +  W|  +1.  (4.7) 

An  approximate  joint  distribution  of  the  transformed  observations  is 
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P  {Zi  £dzi,Z2  ^dz2,-^*,Zj+i  edzj+i} 

=  P  {Zi  €  dZi}P{Z2  €  dZ2\Zi  =  2^}^  —  +;[  €  dZj^i\Zi  =  2]^, ...  ,Zj  =  27} 

=  P{Yi  =  dyi}P{Z2  €d22|Yi  =  i/i}*...  xPjZj+i  £rf27’+i|Yi  =yi,...,Yj  =  yj} 


1  2 
-  //{(«/ W) 


(4.8) 


where 

ftia,W)  =  [a^Ct  +  x^^iXt  +2apiXt+i^^  +  W][k'{amt  +  Xi^lbi)f  +1 
Hence  an  approximate  In-likelihood  function  is  (up  to  addition  of  constants). 


/(a,W)  =  £  -^\nftia,W)  -  ^(z^+i  -  A:(aw^  +  Xf+jb^)) Vf(a,W) 
f-l  ^  ^ 


(4.9) 


Note  that 


k{x)  =  [4e^  +  l|  ; 


Ir  V  .-05  ^ 

k'{x)  =  -[4e^  +  ij  = 


2e^ 


k"ix) 


2e^[2e^+l\ 
+  ‘ 


(4.10) 


Further, 
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(4.11) 


-^ft{a,]N)  =  [k'ianit  + 

-^ft{a,W)  =  [2aQ  +2xt^iPtJqJ^][k'{amt  +Xf+i6,)]^ 


+[a^Ct  +  Xt^itt  +  2aptXt+ijQV^  +  w|2{fc'(amf 
^k"  {am^  +  Xt+ibf)mt.  (4.12) 

Differentiating  (4.9)  with  respect  to  a  and  W  we  obtain 


dl 


^  ^  oa 


—  =  y  - 

da  2  ftia,w) 


1  - 


[zt^l  •  kjamt+xt^ibt)] 

ft{a,W) 


(4.13) 


+(z,+i  -  k(amt  +Xt+ibt))k'{amt  +  Xf+i&()mf^(a,W) 

a 


T 

=  y  -  - - - —ff{a,W)  +  — [2*+i  -  k(amt 

dW  2f^{a,}N)d^N;*  >  >  ^  ^  ^ 


T 

■■  y  -  - - - - —ftia,W) 

2ftia,W)dW 


^ - ftia,W) 

^^^1^)1  (4.14) 

Ma,wf 

2 " 

[zt-^i  -  kjanit  +Xt+ibt)] 
fticc,W) 


\d^l] 

T 

=  y 

1 

A/,(a,W) 

5a 

_da^ 

u 

t=l 

2 

- 

t2 


[k'{amt  +  Xt+ibt)mt]' 


ftia,W) 


(4.15) 


d^l] 

^  1 
=  -y  - 

— /f(a,W) 

aw 

.aV. 

L  2 

f=l^ 

ft{a,W) 

(4.16) 
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d^l 

^  1 
=  -V  — 

dW  da 

dadW_ 

2^  2 

(4.17) 


To  avoid  difficulties  with  a  Newton  procedure  involving  the  restricted 
range  of  W  >  0,  we  will  reparameterize;  ]N  =  e^.  In  this  case 


_a[  ^  dl 
dy  dW 


1  1 

CM 

<0  CO 

1 _ 1 

=  E 

r  d^i  1 
.dW.^J 

’  d^l  ' 

=  p 

d^l 

dady 

Hi 

dWda 

and  W  is  replaced  by  e  in  all  the  expressions  (4.12)  -  (4.17). 

A  Kalman  procedure  with  this  method  of  estimating  a  and  W  is  similar  to 
that  of  Section  3.2  except  that  Step  1  is  replaced  by  the  numerical  solution  of 
the  system  of  equations. 


(4.18) 


(4.19) 


by  a  Newton  procedure  using  Fisher's  scoring  method.  Set  ]N  =  e'  where  y  is 
the  solution.  Occasionally,  there  are  numerical  problems  with  the  Newton 
procedure.  In  these  cases  search  is  used  to  find  estimates  of  a  and  W. 

A  summary  of  the  Newton  (with  default  search)  procedure  is  as  follows. 

The  procedure  starts  with  the  previously  estimated  a  and  ycalled  Uq  and  Yq. 

0.  SetB=0. 
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1. 


Compute  \f{dl,da)  ,  \iidl,dY)  ,  E  \h\hc\[{\i{d^l,da^)),  E 


dy 

a^/'' 


dy 


'  ^  [aaiy]  “o  ^  ^ 


<  0.0001  go  to  2',  otherwise  go  to  2. 


2.  Solve  the  system  of  linear  equations  (in  terms  a  and  f) 


0  =  — +  E 
da 


’^1 

da^ 


{a  ■  «(,)  +  £ 


dh 


dady 


(r  -  ro) 


0  =  — +  E 

ay 


a^/ 

aaay 


(a  -ao)  +  E 


dy^ 


(r  -  Yo) 


for  a„  and  y„.  Go  to  3. 

2'.  Set  y^  =  yQ.  If  ^  is  of  the  same  sign  for  both  the  current  and  previous 
value  of  a,  set 


/r 

a„.ao-—/E 


da^ 


01 

and  go  to  3.  If  ^  is  of  different  signs  for  the  current  and  previous 
values  for  a,  then  a  golden  section  search  is  used  to  solve  the  equation 

da 


set  a„  equal  to  the  first  value  of  a  for  which 


dl 


da 


<  0.1;  if  the  number  of 


iterations  for  this  one-dimensional  search  is  greater  than  50,  go  to  4. 


3. 

3'. 


If  max 
values. 


dl 


da 


dl 


da 


(Yn) 


If  I  I  <5  and  y^  <  5,  go 


<0.1  stop  and  take  a^  and  y^  as  the  estimated 

to  5. 
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3".  If  1  a„  I  >  5  or  >  5,  let  B  =  B  +1.  If  B  =  1,  go  to  4.  If  B  =  2,  go  to  4'.  If 
B  >  3,  to  to  4". 

4.  Do  a  two-dimensional  search  of  the  In-likelihood  function  (4.9) 
parameterized  in  terms  of  a  and  W  for  a  maximizing  value.  The  grid 
for  a  is  [-4,4]  in  steps  of  0.2;  the  grid  for  VV  is  [0,8]  in  steps  of  0.2.  If 


max 


a/ 


da 


iL 

ay 


(InVV) 


<0.1 


return  and  y^  =  InW  as  the  estimated  values;  otherwise  set  Uq  = 
and  Yq  =  InVV  and  return  to  1. 


4'.  Same  as  4  except  the  grid  for  a  is  [-8,8]  in  steps  of  0.1  and  the  grid  for  W 
is  [0.1,8]  in  steps  of  0.1  with  the  additional  points  0.000001,  0.00001, 
0.0001,  0.001,  0.01. 


4”.  y  is  set  equal  to  its  previous  value  and  a  golden  section  search  for  the 
zero  of  ^  IS  done  as  in  2  . 


5.  Set  cc^  =  and  y^  =  y„  and  return  to  1  if  the  number  of  iterations  is  less 

than  50.  If  the  number  of  iterations  is  greater  than  50  set  B  =  B  +  1.  If 
B  =  1,  go  to  4.  If  B  >  6,  return  and  y„  as  the  estimated  values.  If 

1  <  B  <  5,  go  to  6. 


dl 

dl  . 

dl ,  . 

dl.  . 

6.  Compute 

and 

.  n  1 

<  1 

otherwise  go  to  6". 


6'.  Fix  =  Oq  and  do  a  search  over  the  interval  [-4B,  4B]  for  that  yfor 
which 


0  = 


ay 


dl 


<  0.1;  if  the  number  of 


Set  y^  equal  to  the  first  value  of  y  for  which 
iterations  is  greater  than  50,  go  to  4. 

Fix  y^  =  yQ  and  do  a  search  over  the  interval  [-4B,  4B]  for  that  a  for 

which  0  =  ^  /*  set  equal  to  the  first  value  of  a  for  which  |  ^  <  0.1; 

if  the  number  is  greater  than  50,  go  to  4. 
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The  following  are  empirical  observations.  The  procedure  involves  a  great 
deal  of  computing.  For  small  times  t,  the  estimate  of  W  can  be  very  close  to 
zero;  a  lower  bound  of  0.1  was  placed  on  the  value  of  VV  that  is  input  into  the 
Kalman  procedure.  Further,  a  lower  bound  of  -2  was  also  put  on  The 


stopping  criteria  of  max 


a/ 


da 


dl 


dy 


<  0.1  seems  to  be  adequate;  a  smaller 


tolerance  can  result  in  the  procedure  becoming  unstable  and  greatly 
increasing  the  computational  effort. 


5.  ESTIMATION  OF  a  AND  W  FOR  THE  POISSON/NORMAL  MODEL: 

LOGARITHMIC  TRANSFORMATION 

Results  of  Lambert  (1989)  suggest  that  when  count  data  arise  from  a 
mixture  of  Poisson  distributions,  then  a  smaller  power  transformation  of  the 
counts  than  the  square  root  is  needed  to  stabilize  the  variance  of  the  count 
data.  In  particular,  if  the  variability  of  the  mixture  is  largish,  then  a  log 
transformation  may  stabilize  the  variance.  In  this  section  simple  procedures 
for  estimating  a  and  W  based  on  a  log  transformation  of  the  count  data  are 
described. 

Suppose  Yj  and  p,  are  as  in  (1.1)  and  (1.2).  Given  p,  the  first  two  terms 
of  a  Taylor  expansion  yield 

lnY,^l  »  ln//,^l  +  +  pj  +1  +  [htexp{Xt^i^  +  -  /i/exp{x,^iP  +  p,+i}] 

(5.1) 


Let 


Z 


/+! 


ln//,+i  - 


(5.2) 


Using  the  first  term  of  the  Taylor  expansion  (5.1) 
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Zf+1  “^f+1  =  +“,+1  »  am,  +©(+1. 


(5.3) 


Hence,  given  [Yq,  Y^,  ...,  Yj],  the  distribution  of  has  approximate  mean  am, 
and  variance  W.  We  will  approximate  the  distribution  of  with  a  normal 
distribution  having  mean  am,  and  variance  W.  A  In-likelihood  function 
under  this  approximation  is  (up  to  addition  of  constants) 


Thus, 


T-l 

/(a,W)  =  Y,  - 

t=i 


-InW  --(2,+i  -  am,)^— ; 
2  2  ^  W 


T - 1  ^ 

—  =  y  (2,+i  -am,)m, —  =  0; 

aa  ,t'i  ^  '  W 


—  =  Y*-— + 

aw  2W  2'  '  ^ 


-  am,) 


2 


=  0. 


T-l 

^LiT)=^n - 

Zm,' 

(-1 


WLiT)  = 


1 

T  -  1 


T-l 


I 

^  =  1 


(2f+l 


dLiT)mtf 


(5.4) 

(5.5) 

(5.6) 


(5.7) 


(5.8) 


are  the  resulting  estimates  of  a  and  W. 

Simulation  experiments  suggest  that  di  and  tend  to  have  a  large 
positive  bias. 

Another  possibility  is  to  estimate  a  by 

jT-l  T-l 

«m(7')  =  7  Z  /  Z 

^  (=1  f=l 


as  in  (3.7)  of  Section  3  and  W  by 
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1 


(5.10) 


Wm,L(T) 


T  -1 


I  (Zf+i  -  dr„iT)mtf . 
t-i 


These  estimates  are  easy  to  compute.  Simulation  experiments  suggest  that 
this  estimate  of  tends  to  have  a  negative  bias  and  VV^  ^  has  a  positive  bias 

A. 

but  not  as  large  as  that  of 


6.  RESULTS  FROM  SIMULATION  EXPERIMENTS 

In  this  section  we  report  results  of  simulation  experiments  concerning  the 
behavior  of  the  approximate  Kalman  filter  of  Section  2  under  various 
procedures  for  estimating  a  and  W. 

Each  replication  of  the  simulation  consists  of  generating  a  Poisson  time 
series  {y^;  f  =  0, 1,  2, T}  using  equations  (1.1)  -  (1.2).  The  random  numbers 
were  generated  using  LLRANDOM  II,  cf.  Lewis  et  al.  (1981).  If  a  <  1,  then  Pg 
is  drawn  from  the  stationary  distribution  of  f  >  0}.  The  Kalman  procedure 
with  each  of  the  procedures  for  estimating  a  and  W  is  then  applied  to 
{x/y  f  =  0, 1,  2, ...,  T}.  For  each  time  t,  the  estimate 

At  =  exp{fttexp{pt  +  (6.1) 


is  computed  and  the  mean  square  prediction  error 


1^-1  .2 
^  ^  f=l 


(6.2) 


is  calculated.  In  addition  the  average  estimated  values  of  a  and  IV  and  j3  are 
computed: 


Moc)  =  7!  «( 


t»l 


(6.3) 
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(6.4) 


^  f  =  l 
1  ^ 

miP)  =  -Y,h  (6.5) 

where  df  an(i  VV^  are  the  estimates  of  a  and  W  obtained  after  the  observation  at 
time  t  and  b^  is  as  in  Section  1. 

The  simulation  is  replicated  N  times.  If  ej,m^{a),  m^{W),  m0)  represent 
the  summary  statistics  from  the  t  replication,  then  the  mean  of  the  summary 
statistics  are  computed  for  the  N  replications;  i.e. 

1  ^ 

e=  —  Yei;  (6.5) 

^  N 

Vvarle]  =  rr— rX 

^  -i/-! 

1  ^ 

a  =  —  m,(a);  (6.7) 

^  /=l 

-  1  ^ 

W  =  niiiW);  (6.8) 

-  1  ^ 

(6.8) 

^  j-i 

These  latter  statistics  are  then  compared  for  the  different  procedures  of 
estimating  a  and  W. 

Results  for  the  following  procedures  for  estimating  a  and  W  are  reported. 
1.  DD:  a  and  W  are  both  known  and  are  not  estimated. 
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2.  FT:  Both  a  and  W  are  estimated  by  a  two-dimensional  search  of  the  In- 
likelihood  based  on  the  Freeman-Tukey  transformation  (4.9).  The  grid 
for  W  is  [0.1,3]  in  increments  of  0.1;  the  grid  for  a  is  [-3,3]  in  increments 
of  0.1. 

3.  M/FT:  a  is  estimated  by  the  moment  estimator  using  (3.7);W  is 
estimated  by  searching  the  one-dimensional  likelihood  based  on  the 
Freeman-Tukey  transformation  with  d  set  equal  to  (3.7);  the  grid  is 
[0.1,3]  with  steps  of  size  0.1. 

4.  M/LN:  a  is  estimated  using  the  moment  estimate  of  (3.7).  W  is 
estimated  using  the  moment  estimator  (5.10)  based  on  the  logarithmic 
transformation. 

5.  LN/LN:  a  and  W  are  estimated  using  the  moment  estimators  (5.9)  and 
(5.10)  based  on  the  logarithmic  transformation. 

All  the  procedures  have  a  lower  bound  of  -2  on  a  lower  bound  of  0.1 
on  VVj,  and  an  upper  bound  of  1  on  the  absolute  value  of  d^ 

For  the  results  of  the  simulation  experiments  reported  in  Table  1,  a  =  0.5, 
W  =  0.25,  Xf  =  l,h^  =  1,  P  =  0.5.  The  initial  values  of  the  estimates  VVq  =  0.25,  &q  = 

0.5,  fiQ  =  0;  for  the  Kalman  Cq  =  1,  tg  =  1,  /5q  =  0,  ^  =  0.  The  two-dimensional 

search  for  the  estimates  of  a  and  W  in  FT  takes  a  large  amount  of  computer 
time.  As  a  result,  the  number  of  replications  for  the  experiments  reported  in 
Table  1  is  small.  However,  the  results  of  Table  1  suggest  that  the  two  most 
promising  procedures  are  M/LN  and  M/FT. 

In  Table  2,  the  number  of  replications  is  100.  The  procedures  M/LN  and 
M/FT  only  are  compared.  In  Table  2,  the  initial  values  of  VVg  =  0.25,  dg  =  0.5, 

/ig  =  0;  for  the  Kalman  Cg  =  1,  tg  =  1,  ^g  =  0,  ^  =  0.  . 

In  Table  3,  (xj  =  {0.25,  0.5, 1, 0.25,  0.5, 1, ...},  /3  =  0.5,  W  =  0.25,  a  =  0.5.  The 
initial  values  of  VVg  =  0.25,  dg  =  0.5,  /ig  =  0;  for  the  Kalman  Cg  =  1,  tg  =  1,  /5g  =  0, 

^  =  0. 

In  Table  4,  the  parameters  are  the  same  as  in  Table  3  except  W  =  1.  The 
initial  values  of  VVg  =  0.25,  dg  =  0.5,  fiQ  =  0;  for  the  Kalman  Cg  =  1,  fg  =  1,  /5g  =  0, 

^  =  0. 
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Comparison  of  Tables  2  and  3  suggests,  not  surprisingly,  that  the  mean 
square  prediction  error  is  smaller  when  there  is  a  variable  covariate  {Xj}. 

Comparison  of  Tables  3  and  4  suggests,  not  surprisingly,  that  a  larger 
value  of  W  results  in  a  larger  mean  square  prediction  error. 
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TABLE  1. 


SERIES  LENGTH  =  5;  20  REPLICATIONS 


METHOD: 

DD 

FT 

M/FT 

M/LN 

LN/LN 

Mean  MSE 

4.34 

4.91 

4.77 

4.19 

4.96 

St.  Dev.  MSE 

3.21 

3.82 

3.64 

3.17 

3.76 

Max  MSE 

11.0 

13.2 

11.7 

12.5 

13.5 

Mean  d 

- 

0.10 

0.11 

0.15 

0.26 

Mean  VV 

- 

0.48 

0.67 

0.98 

0.77 

SERIES  LENGTH  =  10; 

20  REPLICATIONS 

METHOD: 

DD 

FT 

M/FT 

M/LN 

LN/LN 

Mean  MSE 

2.94 

3.15 

3.08 

2.77 

3.36 

St.  Dev.  MSE 

2.27 

2.42 

2.28 

1.93 

2.54 

Max  MSE 

10:4 

11.1 

10.1 

8.27 

11.2 

Mean  d 

- 

0.12 

0.08 

0.09 

0.19 

Mean  VV 

- 

0.23 

0.32 

0.78 

0.69 

SERIES  LENGTH  =  20; 

40  REPLICATIONS 

METHOD: 

DD 

FT 

M/FT 

M/LN 

LN/LN 

Mean  MSE 

3.38 

3.90 

3.72 

3.29 

3.97 

St.  Dev.  MSE 

2.00 

2.81 

2.43 

2.16 

2.62 

Max  MSE 

10.7 

13.5 

12.1 

10.3 

11.2 

Mean  d 

- 

0.33 

0.15 

0.17 

0.31 

Mean  VV 

_ 

0.22 

0.31 

0.79 

0.73 
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TABLE  2.  x^sl,a  =  0.5,  W  =  0.25 


SERIES  LENGTH  =  5;  100  REPLICATIONS 


METHOD: 

DD 

M/FT 

M/LN 

Mean  MSE 

4.21 

4.75 

4.32 

St.  Dev.  MSE 

4.68 

5.81 

5.40 

Max  MSE 

28.7 

38.0 

34.1 

Mean  & 

- 

0.08 

0.09 

Mean  W 

- 

0.55 

0.93 

Mean  ^ 

0.34 

0.25 

0.23 

SERIES  LENGTH  =  10; 

100  REPLICATIONS 

METHOD: 

DD 

M/FT 

M/LN 

Mean  MSE 

3.98 

4.41 

4.14 

St.  Dev.  MSE 

3.27 

3.87 

3.58 

Max  MSE 

16.4 

19.0 

17.3 

Mean  & 

- 

0.07 

0.10 

Mean  W 

0.45 

0.88 

Mean  p 

0.42 

0.36 

0.31 

SERIES  LENGTH  =  20; 

100  REPLICATIONS 

METHOD: 

DD 

M/FT 

M/LN 

Mean  MSE 

3.96 

4.35 

4.18 

St.  Dev.  MSE 

2.59 

3.30 

3.06 

Max  MSE 

20.3 

26.5 

23.6 

Mean  & 

- 

0.10 

0.13 

Mean  IV 

- 

0.41 

0.82 

Mean  ^ 

0.48 

0.43 

0.37 
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TABLE  3.  VARIABLE  {x^],  a  =  0.5,  W  =  0.25 


SERIES  LENGTH  =  5;  100  REPLICATIONS 


METHOD: 

DD 

M/FT 

M/LN 

Mean  MSE 

3.69 

3.89 

3.70 

St.  Dev.  MSE 

3.83 

4.09 

4.00 

Max  MSE 

15.6 

18.4 

19.9 

Mean  & 

- 

0.13 

0.12 

Mean  VV 

- 

0.52 

0.80 

Mean  ^ 

0.32 

0.21 

0.18 

SERIES  LENGTH  =  20; 

100  REPLICATIONS 

METHOD: 

DD 

M/FT 

M/LN 

Mean  MSE 

3.11 

3.32 

3.21 

St.  Dev.  MSE 

2.48 

2.80 

2.68 

Max  MSE 

20.4 

22.7 

20.2 

Mean  & 

- 

0.12 

0.10 

Mean  W 

- 

0.35 

0.78 

Mean  p 

0.41 

0.36 

0.28 
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TABLE  4.  VARIABLE  {Xf},  a  =  0.5,  W  =  1 


SERIES  LENGTH  =  5;  100  REPLICATIONS 


METHOD: 

DD 

M/FT 

M/LN 

Mean  MSE 

19.21 

20.70 

19.88 

St.  Dev.  MSE 

52.2 

53.0 

52.1 

Max  MSE 

415.1 

390.1 

393.5 

Mean  & 

- 

0.19 

0.18 

Mean  W 

- 

0.93 

1.05 

Mean  ^ 

0.31 

0.14 

0.12 

SERIES  LENGTH  =  20; 

100  REPLICATIONS 

METHOD; 

DD 

M/FT 

M/LN 

Mean  MSE 

18.93 

21.6 

20.4 

St.  Dev.  MSE 

28.6 

31.5 

30.9 

Max  MSE 

181.7 

192.1 

188.5 

Mean  & 

- 

0.20 

0.19 

Mean  VV 

- 

1.2 

1.1 

Mean  ^ 

0.55 

0.32 

0.31 
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Of  the  two  procedures,  M/LN  tends  to  have  the  smaller  mean  MSE. 
M/LN  tends  to  have  a  positive  bias  estimating  W  and  negative  biases 
estimating  6.  and  p.  The  procedure  M/FT  takes  more  computational  effort, 
tends  to  underestimate  W,  a  and  P  and  produces  slightly  higher  mean  MSE. 

Tables  5-7  report  results  comparing  procedures  estimating  a  and  W  with 
the  procedure  NET  in  which  both  oc  and  W  are  estimated  using  the  Newton 
procedure  of  Section  4.  For  all  the  procedures  reported  in  these  tables,  there 
is  no  bounding  on  &,  In  Table  5,  x^=  1,  h^•  1,  P  =  0.5,  a  =  0.5,  W  =  0.25.  In 
Table  6,  x^  =  0.25, 0.5,  1, 0.25, 0.5,  1, ...,  and  the  other  parameters  are  as  before. 
For  Tables  5  and  6  the  initial  values  of  VVq  =  0.25,  =  0.5,  pQ  =  0,  Cq  =  1,  fQ  =  1, 

pQ  =  0.  The  values  in  parentheses  are  estimates  of  standard  deviations;  e.g. 

r  1  ^  _2i2 

For  Table  7,  {xj  is  as  in  Table  6,  a  =  0.5,  W  =  0.25,  Pq  =  0.5,  and  the  initial 
estimates  are  &q  =  0,Wq  =  1,  Cq  =  1,  =  1,  pQ  =  0,  /2q  =  0. 

The  Newton  procedure  of  Section  4  does  not  take  as  much  time  as  the 
two-dimensional  search  used  in  Tables  1-2;  however,  it  is  still  a  much  larger 
computational  effort  than  the  other  two  procedures.  Once  again,  based  on  the 
mean  MSE  and  computational  effort  the  procedure  M/LN  appears  to  be  the 
most  attractive.  The  procedure  M/LN  once  again  appears  to  overestimate  W 
and  underestimate  a  and  p. 
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TABLE  5.  X,  s  1;  a  =  0.5,  W  =  0.25,  =  0.5 

SERIES  LENGTH  =  5,  25  REPLICATIONS 


METHOD: 

DD 

NET 

M/ET 

M/LN 

Mean  MSE 

3.94 

4.44 

4.32 

3.89 

Std  Dev.  MSE 

3.08 

3.66 

3.49 

3.10 

Max  MSE  . 

11.0 

13.2 

11.7 

12.5 

Mean  & 

- 

0.06  (0.37) 

0.08  (0.22) 

0.10  (0.23) 

Mean  VV 

- 

0.46  (0.47) 

0.60  (0.60) 

0.93  (0.45) 

Mean  ^ 

0.28 

0.18 

0.21 

0.18 

SERIES  LENGTH  =  20,  25  REPLICATIONS 

METHOD: 

DD 

NET 

M/ET 

M/LN 

Mean  MSE 

3.11 

3.40 

3.32 

3.14 

Std  Dev.  MSE 

1.83 

2.12 

2.11 

2.03 

Max  MSE 

8.50 

9.29 

9.60 

8.78 

Mean  d 

- 

0.22  (0.39) 

0.10  (0.30) 

0.12  (0.29) 

Mean  VV 

1 

0.29  (0.24) 

0.32  (0.26) 

0.80  (0.23) 

Mean  P 

0.38 

0.23 

0.31 

0.23 
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TABLE  6.  VARIABLE  {ar^} ;  a  =  0.5,  W  =  0.25,  p  =  0.5 


SERIES  LENGTH  =  20,  25  REPLICATIONS 


METHOD: 

DD 

NET 

M/FT 

M/LN 

Mean  MSE 

2.63 

2.77 

2.81 

2.63 

Std  Dev.  MSE 

1.78 

1.95 

2.06 

1.79 

Max  MSE 

9.24 

10.10 

10.86 

9.12 

Mean  & 

- 

0.21  (0.38) 

0.06  (0.32) 

0.06  (0.34) 

A 

Mean  W 

- 

0.20  (0.14) 

0.25  (0.20) 

0.68  (0.18) 

Mean  P 

0.33 

0.21 

0.27 

0.18 
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TABLE  7.  VARIABLE  {x^} ;  a  =  0.5,  W  =  0.25,  =  0.5 

6Cq  =  0,  Wq  =  1,  C  q  =  1,  •?()  =  1,  /5q  =  0,  /][q  =  0 

SERIES  LENGTH  =  20,  100  REPLICATIONS 


METHOD: 

DD 

NET 

M/FT 

M/LN 

Mean  MSE 

2.95 

3.25 

3.21 

3.04 

Std  Dev.  MSE 

2.21 

2.70 

2.67 

2.28 

Max  MSE 

13.2 

16.5 

16.5 

13.8 

Mean  & 

- 

0.12  (0.34) 

0.11  (0.28) 

0.09  (0.29) 

Mean  VV 

- 

0.34  (0.44) 

0.36  (0.32) 

0.73  (0.19) 

Mean  ^ 

0.36 

0.25 

0.28 

0.21 
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7.  SUMMARY  AND  CO  NICLUSIONS 

In  this  paper  we  have  investigated  approximate  methods  for  estimation 
and  prediction  for  the  hierarchical  Poisson/normal  time  series  model  given 
in  (1.1)-(1.2).  For  given  values  of  the  random  walk  parameters,  a  and  W,  the 
joint  distribution  of  (P,  is  approximated  by  a  bivariate  normal  distribution 
using  the  Laplace  method.  Various  heuristic  methods  for  estimating  a  and  W 
are  presented.  Based  on  simulation  results  and  ease  of  computation,  it  is 
recommended  that  a  be  estimated  by  (3.7) 

T  T 

d(T)  = 

f=i  f'l 

and  W  be  estimated  by  (5.10) 

1 

W{T)=  —  '£{zt,i-aiT)mtr 

^  ^  f-1 

where  is  given  by  (5.2). 

These  estimates  of  a  and  W  along  with  the  approximate  Kalman 

procedure  using  the  Laplace  method  provide  a  computationally  easy 

procedure  for  prediction  of  the  time  series. 

Figures  4-7  show  results  of  a  simulation  of  model  (1.1)-(1.2)  for 

t  =  1,  ...,  100.  In  the  example  fi  =  0.5,  {xj  takes  the  values  {0.25,  0.5, 1}  over  and 

over,  and  =  1.  {cOj}  are  iid  normal  with  mean  0,  variance  0.25,  and  a  =  0.5. 

The  simulation  starts  with  drawn  from  the  stationary  distribution  of  (pj,  a 

2  1 

normal  distribution  with  mean  0  and  variance  W/(l-a  )  =  3-  The  values  of  a 
and  W  are  estimated  using  (3.7)  and  (5.10).  Initial  values  of  a  and  W  are 
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ttg  =  0.5  and  W  =  1,  Cq  =  1,  Tq  =  1,  Pq  =  0,  j3Q=  0,  Pq  =  0.  There  are  no  bounds  on 
/if,  d  and  VV  . 

Figure  4  displays  the  count  data.  Also  displayed  are  the  true  Af  (circles), 
the  estimated  Xf  when  both  a  and  W  are  known  (dashed  line)  and  the 
estimated  Xf  when  both  a  and  W  are  estimated  (solid  line).  The  A,  when  both  a 
and  W  are  estimated  is  perhaps  a  little  more  responsive  to  changes  in  the  data 
than  when  a  and  W  are  known.  The  difference  between  the  two  estimates  of  Xf 
is  greatest  for  small  times  t. 

Figure  5  presents  plots  of  the  estimated  y/W  and  a.  Note  that  W  has  a 
positive  bias  which  will  make  the  Kalman  procedure  more  responsive  to  the 
count  data.  The  estimated  values  of  a  appear  to  be  reasonable. 

Figure  6  presents  estimates  of  /3  and  its  estimated  standard  deviation 
both  for  the  Kalman  with  a  and  W  known  (dotted  line)  and  with  a  and  W 
unknown  (solid  line).  The  estimates  of  Tj  generally  decrease  as  t  increases 
reflecting  the  model  assumption  that  j3  is  constant.  The  estimates  of  ^  appear 
to  be  reasonable.  The  estimates  of  p  with  a  and  W  also  estimated  tend  to  be 
smaller  than  those  for  which  a  and  W  are  known;  this  behavior  may  be  due  to 
the  fact  that  estimation  of  a  and  W  is  accounting  for  some  of  the  variability  of 
the  data  that  would  otherwise  be  accounted  for  by  p. 

Figure  7  displays  the  true  Pf  (circles),  the  estimated  /f,  =  m  f  with 
parameters  a  and  W  known  (dotted  line)  and  the  estimated  /if  with  parameters 
a  and  W  estimated  (solid  line).  The  count  data  are  also  displayed.  The 
estimates  of  /if  with  a  and  W  estimated  are  more  variable  than  those  when  a 
and  W  are  known  reflecting  the  greater  responsiveness  of  the  Kalman  to  the 
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data  when  the  estimated  W  has  a  position  bias.  Once  again  the  estimates  of 
appear  to  be  reasonable. 
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